function t54
global g k m1 m2 h nt x
g=10; m1=1; m2=1; k=100; h=2;
opts=odeset('events',@myevents);
tstart=0;
tfinal=10;
tspan=tstart:0.005:tfinal;
y0=[h,0,0,0];
tout=[];
xout=[];
t=0;
while t(end)<tfinal-0.1%for i=1:100
    
    [t,x,te,xe,ie]=ode45(@f,tspan,y0,opts);
    nt=length(t);
    tout=[tout;t];
    xout=[xout;x];
    tstart=t(nt);
    tspan=tstart:0.005:tfinal;
    Y=fsolve(@qg,[0 0],optimset('Display','off'));
    y0=[x(nt,1),Y(1),x(nt,3),Y(2)];
end

figure
axis equal;
axis([-1 1 -3 3])
lx=[-0.2,0.2];
ly=[0,0];
horiline=line(lx,ly,'color','g','linewidth',2,'erasemode','xor');
theta=0:0.1:2*pi;
xb=0.1*cos(theta);
yc=0;
yb=yc+0.1*sin(theta)+0.1;
ball=line(xb,yb,'color','b','linewidth',2,'erasemode','xor');
yv1=(0:0.1:20*pi)/(20*pi);
yv=yv1*3-3;
xv=0.1*sin(0:0.1:20*pi);
vertiline=line(xv,yv,'color','r','erasemode','xor');
ll=length(tout);
for i=1:ll
ly=[xout(i,3),xout(i,3)];
yc=xout(i,1);
yb=yc+0.1*sin(theta)+0.1;
yv=yv1*(3+xout(i,3))-3;
set(horiline,'XData',lx,'YData',ly);
set(ball,'XData',xb,'YData',yb);
set(vertiline,'XData',xv,'YData',yv);
drawnow;
end
plot(tout,xout(:,1),tout,xout(:,3));
function Y=qg(X)
    Y(1)=X(1)-X(2)+x(nt,2)-x(nt,4);
    Y(2)=m1*X(1)+m2*X(2)-m1*x(nt,2)-m2*x(nt,4);
end
    function [value,isterminal,direction] = myevents(t,x)
% Locate the time when the ball droping on the slice
% and stop integration.
value = x(1)-x(3); % detect height = x(1)-x(3)
isterminal = 1; % stop the integration
direction = -1; % negative direction
    end
        function F=f(t,x)
            
        F=[x(2);
        -g;
        x(4);
        -k/m2*x(3)];
        end
end